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We show how the sign problem occurring in dynamical simulations of random matrices at nonzero 
chemical potential can be avoided by judiciously combining matrices into subsets. For each subset 
the sum of fermionic determinants is real and positive such that importance sampling can be used 
in Monte Carlo simulations. The number of matrices per subset is proportional to the matrix 
dimension. We measure the chiral condensate and observe that the statistical error is independent 
of the chemical potential and grows linearly with the matrix dimension, which contrasts strongly 
with its exponential growth in reweighting methods. 
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Dynamical simulations of QCD at nonzero chemical 
potential are hampered by the sign problem (see ref. [Ij 
for a review). The problem occurs when the probabilistic 
weight used in Monte Carlo methods becomes complex, 
thus precluding the use of standard importance sampling 
techniques. Methods to study QCD at small chemical po- 
tential, where the sign problem is mild, include reweight- 
ing methods, Taylor expansions and analytic continua- 
tion from imaginary chemical potential [1] . QCD at zero 
and nonzero density can also be investigated using ran- 
dom matrix theory because of the universality proper- 
ties of the Dirac operator to leading order in the epsilon 
regime. Based on this equivalence, the severity of the 
sign problem has been studied using the average phase 
of the fermion determinant . As dynamical simula- 
tions of random matrices at nonzero chemical potential 
also suffer from the sign problem they can be used as a 
playground for its study. In this letter we present a sub- 
set method which avoids the sign problem in this case. 

In chiral random matrix theory at nonzero chemical 
potential fj, the Dirac operator for a fermion of mass m 
can be represented by [7| 
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where the configurations = (0i,02) are pairs of com- 
plex random matrices (pi and 02 of dimension (N+iy) x N 
with ly the number of zero modes of D (for m = 0). 
The partition function of the chiral random matrix model 
with Nf dynamical quarks with equal mass m is 

d(bid^2 w;(</)i)u;(02) det^^i?^,™(0), (2) 

with Gaussian weights 

w{<P,) = (iV/^)^(^+'^) exp(-iV tr 4<^,) (3) 

and integration over the real and imaginary parts of all 
matrix components. In the presence of a chemical po- 
tential the fermion determinant becomes complex, the 



sign problem arises and the work needed to make reli- 
able measurements on the statistical ensemble grows ex- 
ponentially with the volume. In reweighting methods this 
exponential increase comes from the need to compute ex- 
ponentially small reweighting factors from a statistical 
sampling of largely canceling contributions [T] . 

Herein we propose a method to avoid the sign problem 
in dynamical simulations of random matrices. The idea 
is to rewrite the partition function ^ as an integral over 
subsets n = : i — l,...,Ns}, each containing Ns 
configurations — {(l)\,(j)2), such that 



(4) 



where by construction the Gaussian weights w{4>\)w{(f)2) 
will be invariant for all 0* in the subset and can thus 
be denoted by W{ft), and we define the fermionic subset 
weight 
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which depends on m and /i. The construction of the 
subsets will be described below, but its crucial property 
is that the sum ai^^mi^) over all the determinants in a 
subset is real and positive, for every fl of the ensemble, 
and can thus be used as a probabilistic weight in a Monte 
Carlo method. The sample average of an observable O 
over a sample of A'mc subsets flk is computed as 
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where (j)'''^ = ((/i^^^l*) G ilk- Methods based on partial 
integrations/summations have also been considered by 
other authors, see, e.g., refs. [SHU]- 

We now sketch how the subsets are constructed. As 
the partition function ^ is real and positive (for < 1) 
the basic idea is to create subsets of matrices for which 
a^^m has the same property. For /i = the Dirac ma- 
trices are anti-Hermitian and their determinants are real 
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and positive. When /i is increased the matrices become 
non-Hermitian, the determinants complex and the sign 
problem sets in. For ii — 1 and m — the ensemble 
becomes maximally non-Hermitian, the average phase 
factor is zero and the sign problem maximal. For this 
limiting case we ideally want to find a partner configura- 
tion ip = (tpi,ip2) for any given (/) of the ensemble, such 
that their fermionic determinants explicitly cancel. This 
is easily achieved if, given some angle 9, we can find a 2p 
for which i3i,o(V') = e^^i^i^ol^)- Using ([T]) the latter can 
be written as 

where 

V'i('^; 6*) = cos6' (^1 -I- sin6' 02 
■i/'2 (0; ^) — cos 9 (j)2 — sin 9 0i 

are orthogonal rotations of {4>i,4'2)i which emerge when 
e^^ Difi{(f)) is explicitly rewritten in an anti-Hermitian 
and a Hermitian part. From ([s]) it follows that 
w{ipi)'w{ip2) = w{4'i)w{4'2) such that ip and have 
the same Gaussian weights in the partition function, for 
any rotation 9. An exact pairwise cancellation of the 
fermionic determinants will happen for any 9 obeying 
^i2NNfe _ rpj^g jjgg^ ^£ canceling determinants can be 
extended from pairs {(pjip) to subsets {'ilj{(j>;9n)}, where 
the 9n are such that the total sum of fermionic determi- 
nants cancels for ^ = 1 and m = 0, i.e., gi2NNf0„ _ q 
This idea is now ported to arbitrary chemical potential 
/i and mass m. In the following we will consider subsets 

= 0„) : 0„ = — A n - 0, . . . ,iV,-l j , (9) 

with -0 = ("01, "02) defined in ([s]) and Ng a positive integer. 
To each configuration = (0i,02) of the random ma- 
trix ensemble corresponds a subset ^2(0), and the set of 
all subsets forms an TVg-fold covering of the random ma- 
trix ensemble. The configurations '0(0i ^) have Gaussian 
weights independent of 9, but different fermionic weights 
det ^^-0^,^(0(0;^)). 

The method presented in this paper is based on the 
following theorem: For any Q constructed according to 
([9]), and for arbitrary /i < 1 and m, the fermionic sub- 
set weight crp „i(f2) defined in ([s]), i.e., the sum of the 
fermionic determinants of the configurations ■0(0^ ^n)j 
is real and positive if Ns > NNf. 

More specifically, for 771 = it can be shown that 

a^,o{^) = {l-^iX^fao,om, (10) 

for any fi. As cro,o is real and positive, a^^ is also real and 
positive for fj, < 1. For /i = 1 the sum of determinants is 
exactly zero. For nonzero mass we can show that cr^.m 
is real and (T^,m(rj) > (1 — fi^)^^f ao,,n{^) for ^ < 1, 



such that the weight cr^^rn is positive. The weights cr^^m 
can thus be used to generate subsets of random matrices 
using a Metropolis algorithm and compute observables 
using ([6]). In practice, Ng will be set to its minimum 
value, i.e., = NNf + 1. 

This theorem can be proven analytically and was thor- 
oughly tested numerically. The proof will be given in a 
forthcoming paper. 

We applied the subset method to compute the chi- 
ral condensate of the chiral random matrix ensemble, 
which is defined as E = tr [T3|. We use the 
Metropolis algorithm to generate subsets O according to 
their weights VF(f2)(T^.,„(ri), where the determinants are 
computed numerically. Successive subsets in the Markov 
chain are generated as follows: randomly choose a con- 
figuration in the current subset, generate a new config- 
uration by making a random step, construct the subset 
corresponding to this new configuration, apply an accept- 
reject step to the newly proposed subset. In each Markov 
chain we generate 100,000 subsets. The chiral condensate 
is computed for each configuration and its average is com- 
puted using ([g]). As usual, successive configurations in 
the Markov chains are correlated and the number of in- 
dependent subsets is smaller by a factor 2r, where r is the 
integrated autocorrelation time. The statistical errors on 
the measurements are determined using the standard er- 
ror formula corrected for the autocorrelations. 

To compare the subset method with standard reweight- 
ing methods, all the simulations were repeated using 
quenched, phase quenched and sign quenched reweighting 
[6]. In those cases, we used {NNf + 1) x 100, 000 random 
matrices in the Markov chains, such that the total num- 
ber of matrices generated in the reweighting methods is 
the same as in the subset method. For the sake of clarity, 
we only show the results of phase quenched reweighting 
in the figures below, as its results are representative for 
the various reweighting methods. 

We performed simulations for = 2, . . . , 34, choosing 
Nf — 1 and m — 0.1/2N, so that the mass is small w.r.t. 
the magnitude of the smallest eigenvalue. In fig. [l] the 
chiral condensate is shown as a function of the chemi- 
cal potential for matrices with A^ = 2,4,8. We compare 
the results obtained using the subset method with those 
from phase quenched reweighting. We also show the ex- 
act results computed in ref. [13]. The data are shown in 
the top row and the corresponding relative statistical er- 
rors in the bottom row. As the matrix size increases the 
reweighting method fails for smaller and smaller /Lt^ due 
to the sign and the overlap problem. As expected, the 
error of the reweighting method grows exponentially, un- 
til the method fails when the set of sampled matrices no 
longer overlaps with the relevant configurations for the 
given value of /i^ and N. This strongly contrasts with 
the results of the subset method which are reliable up to 
much larger values of /i^ and agree with the analytical 
predictions. Interestingly, the accuracy of the measure- 
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FIG. 1. Chiral condensate E (top row) as a function of the chemical potential /i^ for the subset method and the phase-quenched 
reweighting method for = 2, 4, 8. The full line shows the exact analytical result of ref. [IS]- The reweighting method fails for 
ever smaller [P' when A'^ grows (some data points are negative and are left out of the semi-log plots). The corresponding relative 
statistical error e is shown in the bottom row. For the subset method, the error is independent of the chemical potential. The 
error for the reweighting method grows very rapidly and should only be trusted as long as the method works (see top row). 



ments is independent of the chemical potential. 

We also measured the chiral condensate as a function 
of iV for fixed values of /i^ . The A^-dependence of the rel- 
ative statistical error on the chiral condensate is shown 
in fig. [2] (left: subset method, right: phase quenched 
reweighting) for different values of /i^. For a fixed num- 
ber of subsets the error in the subset method (left panel) 
increases approximately as ^pN and is independent of /i 
(the latter also follows from fig. [T]) . If we fix the num- 
ber of matrices, rather than the number of subsets, the 
error will increase with an additional factor (as the 
subset size itself grows with iV -I- 1), such that the over- 
all relative error will grow linearly with N . Turning the 
argument around, to achieve a constant error the num- 
ber of subsets would have to grow proportionally to iV, 
i.e., the total number of matrices should approximately 
grow as iV^. The right plot shows the same quantity for 
phase quenched reweighting (on semi- log scale). We ob- 
serve that the error grows exponentially with N until the 
method completely fails when the error stagnates around 
one and is no longer reliable. Note that for both methods 
the additional cost for the numerical computation of the 
determinants is proportional to . 

For large values of N or /i^ we observed that the subset 
method breaks down at an A^-dependent value fic{N), be- 
cause the finite precision arithmetic limits the accuracy 
of (Jfi^m, which can lead to a violation of its positivity. 
The breakdown occurs when the sum of the determinants 
in a subset becomes of the order of the accuracy of the 
individual determinants. The approximate value /x^ is 



shown as a function of N in fig. [3] for simulations per- 
formed in double precision arithmetics. The data are 
fitted well by /i^ « 1 — exp(— 38/A^^'^), whose functional 
form is inspired by (10). For large N the fitted curve 
goes like The breakdown can easily be detected 

during simulations by monitoring cr^.m and the magni- 
tude of the individual determinants. When the ratio of 
both quantities is of the order of the machine precision 
the method no longer gives sensible results. Usually this 
is accompanied by a numerical violation of the positivity 
of (J^,rn- For simulations at m = this breakdown can 



be avoided altogether by using the analytic formula ( 10 ) 
to compute the sampling weight ct^^q from cto.o- 

Note that the breakdown caused by the finite nu- 
merical accuracy will equally well show up in standard 
reweighting methods, even if an exponentially large ef- 
fort is invested to keep the statistical error under control. 
This accuracy problem is therefore different from the sign 
problem, as the latter is of a pure statistical nature and 
does not depend on the use of finite-precision or exact 
arithmetic. The former is an additional problem of nu- 
merical nature which can be improved upon if necessary 
by using more sophisticated numerical methods. 

What happened to the sign problem in the subset 
method? In standard reweighting methods the large can- 
cellations, inherent to simulations at real chemical poten- 
tial, happen through statistical sampling of the partition 
function. In the subset method these cancellations are 
removed from the statistical sampling procedure and be- 
come confined inside the subsets, which are constructed 
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FIG. 2. Relative error e on the chiral condensate versus matrix size A'^ for various values of chemical potential fi^ = 



0.1,0.2,0.3,0.4,0.5. The left plot shows the results of the subset method. The full curve e{N) oc vTV serves to guide the 
eye. As a comparison, the right plot shows the relative error for phase quenched reweighting on a semi-log plot (the color 
coding for fi^ is the same as in the left plot). The error increases exponentially with A'^, until the reweighting method fails. 




FIG. 3. Chemical potential fic{N) where the subset method 
breaks down due to the finite precision in the computation of 



The full line shows the fit fi^ 



exp(-38/iV' 



in a deterministic way and whose partial sums cr^^m yield 
net real and positive weights. Thus, a fundamental dif- 
ference between both methods is that the number of con- 
figurations in reweighting grows exponentially with the 
volume to maintain the necessary statistical accuracy on 
the average weight factor [T] , whereas only NNf + 1 ma- 
trices per subset are needed in the subset method to com- 
pute the positive subset weights (T^t,m, independently of 
its magnitude. This is how the subset method avoids 
the large statistical cancellations characterizing the sign 
problem. 

From a practical point of view, even though finite- 
precision arithmetic in the numerical simulations even- 
tually leads to a breakdown of the subset method, our 
numerical tests have shown that the method yields a vast 
improvement over the standard reweighting methods for 
the dynamical simulation of random matrices. With the 
new method higher values of /i^ and N can be accessed, 
without any loss of accuracy when increasing the chemi- 



cal potential and with a measurement error growing pro- 
portionally to N when increasing the matrix size. 

The crucial question whether the subset method is also 
applicable to physically relevant systems will be investi- 
gated in future research. 
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by the DFG collaborative research center SFB/TR-55. 
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